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We study the thermalization process in classical Yang-Mills (CYM) field theory starting from noisy 
glasma-like initial conditions by investigating the initial-value sensitivity of trajectories. Kunihiro 
et al. [171 ] linked entropy generation to the Kolmogorov- Sinai entropy, which gives the entropy 
production rate in classical chaotic systems, calculated numerically for CYM fields starting from 
purely random initial field configurations. In contrast, we here study glasma-like initial conditions. 
For small random fluctuations we obtain qualitatively similar results while no entropy increase is 
observed when such fluctuations are absent. We analyze the intermediate time Lyapunov spectrum 
for several time windows and calculate the Kolmogorov- Sinai' entropy. We find a large number of 
positive Lyapunov exponents at the early stages of time evolution. Also for later times their number 
is a sizeable fraction of the total number of degrees of freedom. The spectrum of positive Lyapunov 
exponents at first changes rapidly, but then stabilizes, indicating that the dynamics of the gauge 
fields approaches a steady state. Thus we conclude that also for glasma-like initial conditions a 
significant amount of entropy is produced by classical gluon field dynamics. 



I. INTRODUCTION 



Motivation 



Elucidating the mechanism of entropy production 
and early thermalization leading to quark-gluon plasma 
(QGP) formation is one of the fundamental problems 
posed by the phenomenological analyses of the exper- 
imental data obtained from the Relativistic Heavy-Ion 
Collider (RHIC) at Brookhavcn National Laboratory and 
the Large Hadron Collider (LHC) at CERN. As hydrody- 
namics is only applicable after local thermalization has 
occurred, massive entropy production is a basic ingredi- 
ent for any detailed understanding of the time-evolution 
of heavy-ion collisions. Analyses based on ideal relativis- 
tic hydrodynamic equations suggest that thermalization 
should be achieved in an early stage in order to explain 
the RHIC data. The thermalization time is estimated to 
be To ~ 1 fm/c p], [2[ , which is significantly shorter than 
the equilibration time obtained in perturbative QCD |3[. 
If To were much longer, one would have to expect that 
significant artifacts caused by incomplete thermalization 
could affect many interpretations of the experimental 
data. 

The "early thermalization" of the created matter 
should be explained on the basis of the underlying dy- 
namics of the time-evolution of the matter created in 
the nuclear collision. This remains a major challenge 
of relativistic heavy ion physics. Apart from the time- 
evolution dynamics, peculiar fluctuations in the initial 
conditions are also recognized to be an important ingre- 
dient for understanding the hydrodynamic evolution of 
the quark gluon matter, in particular, the elliptic and 
triangular flows 0t3|- 



B. Classical Yang-Mills Fields 

Classical Yang-Mills (CYM) field theory is a good 
starting point for describing the initial stage dynamics 
of relativistic heavy-ion collisions. At high energies, wee 
partons with small Bjorken x are excited, and the number 
of gluons in the small x region is so high that they may 
be treated as a coherent or classical field (a, Q . A clas- 
sical solution of gluon fields in this regime is known: A 
nucleus at high energy can be regarded as a collection of 
large- a; particles, the valence partons, which act as color 
sources that give rise to strong, transversely polarized 
chromo-magnetic and -electric fields. This configuration 
of the gluon fields is referred to as the color glass conden- 
sate (CGC). When two nuclei collide, their gluon fields 
interact and produce new, longitudinally polarized fields 
that extend between the color sources in the direction of 
the beam axis. This glasma field configuration is used as 
a standard initial condition for the time evolution of the 
matter created in relativistic heavy-ion collisions. 

It is generally expected that some field instabilities 
could cause early thermalization; such instabilities in- 
clude the Weibel instability [T(J, LL1| an d Nielsen- Olesen 
instability |L2| - [l4| . In both cases, the instability gives 
rise to an exponential growth of the amplitude of un- 
stable modes having high energy, which are expected to 
cascade into a large number of modes or particles with 
small energies due to the nonlinear coupling of fields, un- 
til the fields are thermalized. 

There are two facets in thermalization of CYM. The 
first aspect is the isotropization of the energy momen- 
tum tensor. The approximate isotropization is the min- 
imal condition to apply hydrodynamics, and has been 
discussed in the literature [ljl [Tg]. The second is the 
equilibration of the spectrum toward the Bose-Einstein 



and the Fermi-Dirac distribution of gluons and quarks, 
respectively. While the spectrum equilibration provides 
for a more rigorous definition of thcrmalization, it has 
not been discussed seriously so far in CYM. One reason 
for this may be that classical field theories do not gen- 
erally exhibit the correct equilibrium spectrum by them- 
selves. When equipartition of energy is realized in classi- 
cal equilibrium, high momentum components are favored 
compared with the quantum equilibrium, i.e. the Bose 
distribution. In the thermalization scenario described in 
the previous paragraph, we expect the produced particles 
to show quantum equilibrium spectra. Thus the very 
thermalization mechanism is attributed to this conver- 
sion process but it is a difficult and unsolved problem to 
describe non-equilibrium and non-uniform situations in 
terms of nonlinear classical or quantum dynamics. 



C. Entropy Production 

A central task to reach an understanding of early ther- 
malization is to identify the mechanism of entropy pro- 
duction, which has been scarcely touched upon so far in 
the literature. It is to be noted that the high occupation 
probability of the unstable modes makes possible a quasi- 
classical treatment of the thermalization process in the 
initial stage. Furthermore, Kunihiro, Muller, Ohnishi, 
and Schafcr jT7[ noticed that the use of the Husimi-Wehrl 
entropy <Shw, tne Wehrl entropy [18[ defined in terms of 
a smeared Wigner function (the Husimi function) [19| , is 
meaningful in the quasi-classical regime and showed that 
<Shw grows at the rate of the Kolmogorov-Sina'i (KS) 
entropy Sks hi the long-time limit. The fundamental 
origin of entropy growth induced by the Husimi trans- 
formation is that even without an actual measurement 
the quantum mechanical uncertainty principle implies a 
minimum of coarse-graining. In other words: Informa- 
tion which cannot be measured in accordance with the 
uncertainty principle is de facto lost and information loss 
is equivalent to entropy growth. 

Here the KS entropy is defined as a sum of the positive 
Lyapunov exponents Ai, 



'KS 



i,A;>0 



(1) 



The size of a Lyapunov exponent A^ is an indicator for 
the initial- value sensitivity of trajectories defined through 
the equation, \5Xi(t)\ ~ \SXi(to)\ exp[A^(i — t )], where 
SXi represents a small difference of phase-space variables 
at the initial time t = to- A positive Lyapunov exponent 
means that the distance between the two phase space 
points grows exponentially in time. Thus the entropy 
production in classical dynamics is closely related to the 
chaoticity of the system. It is therefore quite pertinent 
to investigate possible entropy production of the clas- 
sical field itself to study the thermalization mechanism 
through its chaotic behavior in the glasma stage. 



In fact, the study of the chaotic properties of the classi- 
cal evolution of Yang-Mills fields has a long history [20] , 
initiated by the observation of chaotic behavior in the 
infrared limit of Yang-Mills theories [21( , which was con- 
firmed later for the compact lattice formalism of classical 
Yang-Mills theory [22| . Some properties of Lyapunov ex- 
ponent and KS entropy in compact lattice gauge theories 
are discussed in [23[ and [2J], respectively. 

In Ref. [25||, the authors analyzed the exponential 
growth of the distance between two trajectories for classi- 
cal Yang-Mills evolution, starting from adjacent generic 
random initial gauge fields in the noncompact {A, E) 
scheme. They calculated the KS entropy and found that 
the KS entropy is positive and finite even after a long 
time. The equilibration time scale r oq was estimated 
to be around 2 fm/c for T = 350 MeV, though with 
rather substantial systematic uncertainties. This result 
obtained by starting from a generic random initial fields 
shows that a significant amount of entropy is produced 
by the dynamical complexity inherent in the CYM equa- 
tions, suggesting that thermalization in heavy-ion colli- 
sions can be at least partly be achieved in the classical 
regime before particle production comes into play as a 
quantum process. 

In Ref. [25[, it was also observed that there are three 
distinct time regimes, namely a kinetic stage for short 
sampling times, an intermediate- and a long-time regime: 
i) The short sampling time is characterized by the lo- 
cal Lyapunov exponents (LLE). ii) The evolution of the 
distance on (long) mixing time scales is described by the 
usual Lyapunov exponents, which we refer to as global 
Lyapunov exponents (GLE). It is to be noted that GLEs 
are equivalent to the original Lyapunov exponents, iii) In 
the intermediate-time period, the nonlinear coupling be- 
tween different field modes is significant but the energy 
remains localized among the primary unstable modes. 
The Lyapunov exponents characterizing the exponential 
growth of the separation of trajectories in this interme- 
diate time period are called intermediate Lyapunov expo- 
nents (ILEs). We emphasize that the ILEs are the most 
relevant Lyapunov exponents for the thermalization of 
the glasma, because they characterize the time-evolution 
of the strongly excited Yang-Mills fields in the early stage 
when the field configuration is still far away from equi- 
librium and a quasi-classical description of the dynamics 
of the Yang-Mills field is appropriate. 



D. Scope of this Work 

In this paper, we extend the analysis done in [25| and 
explore the thcrmalization process starting from more re- 
alistic initial conditions than those adopted in Ref. |25| . 
namely glasma-like initial conditions. We analyze the 
time-evolution of the distance between two trajectories 
starting from two adjacent points in phase space and also 
extract the spectrum of Lyapunov exponents from the 
time evolution of these fields. As was mentioned before, 



we can determine whether and how early entropy produc- 
tion is achieved by examining the number and magnitude 
of the positive Lyapunov exponents. In this study, we use 
two glasma-like initial conditions. One is called "modu- 
lated initial condition" , where the initial color-magnetic 
fields Bi are spatially modulated along the z and x-axes. 
We call the other "constant- A initial condition". Here 
both the gauge potentials Ai and the chromomagnetic 
fields Bi are constant, but non-commuting. It turns out 
that the two initial conditions give similar results for 
entropy production with some minor differences in the 
time evolution once the glasma-like initial conditions are 
taken. 

The paper is organized as follows. In Sec. HH we intro- 
duce the basic ingredients of our simulations including 
the initial conditions. In Sec. IIII1 we show the numeri- 
cal results for calculations starting from the modulated 
initial condition. From the time evolution of the dis- 
tance between two trajectories which are very close to 
each other initially we determine the Lyapunov expo- 
nents. We also show numerical results for simulations 
using the constant- A initial condition in Sec. IIIII Section 
IIVI is devoted to a summary and concluding remarks. 



II. FORMULATION OF THE PROBLEM 

In this section, we first introduce the equations of mo- 
tion and the formulae needed for the analysis of chaotic 
behavior and the entropy production. (See Ref. [25( for 
details.) Then we describe the glasma-like initial condi- 
tions with fluctuations as well as the parameters chosen 
in our simulations. 



A. Classical Yang-Mills equation 

In pure Yang-Mills theory in temporal gauge A$ = 0, 
the Hamiltonian in the noncompact (A, E) scheme takes 
the following form on a cubic spatial lattice 



of motion (EOM) for (Af{x), Ef(x)) 



h =\T, e ^ x ) 2 + \ E W > ( 2 ) 

x,a,i x,a,i,j 

F«(x) = d t A«{x) - djAfix) + £ r bc Al(x)A](x) , 



b.r 



(3) 



where d{ is the central difference operator in the i- 
direction, d t A(x) = {A(x + i)-A(x-i)}/2, and f abc are 
structure constants. In this study, we deal with SU(2) 
gauge theory, i.e., f abc = e abc , where e abc is the Levi- 
Civita tensor defined with e 123 = 1. Note that all quan- 
tities in the equations are dimcnsionless, i.e. scaled with 
appropriate powers of the lattice constant. 

From Eqs. © and (|3|), we get the classical equations 



A?(x)=Ef(x), 

E?{x)=Y,d 3 F^ 



b,c,j 



f abc A b (x)F^x) , 



(4) 
(5) 



which we solve with glasma-like initial conditions to be 
specified later. A fourth-order Runge-Kutta method is 
adopted to solve these EOM. It should be noted that we 
chose the initial condition such that it satisfies Gauss' law 
(EtDi&Wr = E i d i E i -(x)+j: i ^ b r cb A9E ib (x) = 
and check its validity carefully at every time slice. Energy 
conservation is also checked for each time step. 



B. Distance between two trajectories 

A natural indicator for the chaotic behavior of the 
system is the "distance" between two trajectories which 
start from slightly different initial points in phase space. 
With (Af(t, r), Ef(t, r)) being the trajectory starting 
from an initial point (Af(0, r), Ef(Q, r)) in phase space, 
we consider a second trajectory (A'^lt, r), S|°(i, r)) 
starting from an adjacent point. Then we define the dis- 
tance Dee (Dff) between the electric fields (the field 
strengths) by 



D 



EE 



\ x I a,i a,i 



d 



FF 



\ 



E E^(*) 2 -E^( a 



(6) 



(7) 



x I a,i,j 



respectively Here F^{x) is the field strength tensor 
evolved from the initial point (A'°-(t = 0,r), E'°-{t — 
0, r)). These distances are gauge invariant under the 
residual gauge transformation of the EOM, namely, E — > 
D,(x)E and F — > D,(x)F, where Q(5f) is a function 
of gauge transformations in the adjoint representation 
which are independent of time. 

The Lyapunov exponents are extracted from the time- 
dependence of these distances, as is described in the next 
subsection. 



C. Lyapunov exponents 

For the two trajectories (Af(t, x), Ef{t, x)), and 
(Af(t, x), Ef(t, x)) the tangent vector 



6X(t) = (6AZ(t,x),5E°(t,x)) T 
satisfies the following EOMp5J|. 
1 N 



SX(t) 



-H AA {t) 



8X(t) = H(t)5X(t) 



(8) 



(9) 



where we have introduced a matrix defined by 
{H AA (t)) iaxJby = S 2 H/6A*(x,t)8A b Mt). 



(10) 



We call the matrix % Hessian, and the eigenvalues of H 
for each time slice are referred to as the local Lyapunov 
exponents (LLEs) [25| : As it stands, the LLE plays the 
role of a "temporally local" Lyapunov exponent, which 
specifics the departure rate of two trajectories in a short 
time period. 

For a system where stable and unstable modes cou- 
ple with each other as in the present case, an LLE does 
not generally agree with the Lyapunov exponent in a 
long time period. As an adequate way to characterize 
the exponential growth of the fluctuation, the authors in 
Ref. [25J introduced another kind of Lyapunov exponent 
A ILE called the intermediate Lyapunov exponent (ILE), 
which is an "averaged Lyapunov exponent" for an inter- 
mediate time period At; i.e., a time period which is suf- 
ficiently small compared to the thermalization time but 
large enough to sample a significant fraction of phase 
space. 

The explicit definition of a A ILE goes as follows: We 
first note that Eq. © is solved for any time period At 
by 



SX(t + At) = U(t, t + At)SX(t) , 

/ r t+At 



U(t,t + At)=T 



exp 



U(t + t')dt' 



(11) 
(12) 



with 7" denoting the time ordered product. The follow- 
ing Trotter formula for the time-evolution operator U is 
found convenient for a numerical evaluation; 



U(t,t + At)=T Y[ U(t k -i,t k ) 

k=l,N 

-T H [l + H(t k -i)6t] , (13) 



k=l,N 



where St = At/N. Diagonalizing the matrix U, we reach 
the definition of the ILEs 



U D (t,t + At) =diag(e A i 






(14) 



In this study, A ILE 's are calculated by setting t = 0, 

C/ I) (0,Ai)=diag(e Al 1 LEAt ,e^ EA ',...). (15) 



As the thus defined ILEs depend on the time period At, 
we examine the At dependence of the ILE spectrum in 
our later discussion. 

Two comments are in order, here: A Lyapunov ex- 
ponent can be (real) positive, negative, zero or purely 
imaginary. Liouville's theorem tells us that the determi- 
nant of the time evolution matrix U is unity, implying 
that the sum of all positive and negative ILEs is zero. 
The KS entropy is given as a sum of positive Lyapunov 
exponents. The second comment concerns gauge invari- 
ance of the Lyapunov exponents. The discussion in the 
paper is centered on the time evolution of KS entropy 
which is the sum of the positive Lyapunov exponents in 
CYM. However, if the modes related to LLE and ILE 
are gauge non- invariant, the observed chaoticity does not 
necessarily have a physical meaning. In the Appendix we 
show, however, that LLE and ILE are indeed gauge in- 
variant under time- independent gauge transformations in 
the temporal gauge. 

Our goal is to clarify how the entropy grows during the 
time interval when the gauge field configuration is still far 
from equilibrium but has already sampled a significant 
fraction of phase space. Thus the KS entropy of interest 
should be defined as the sum of positive ILEs, 



dS 



— -S KS - ^ A 



ILE 



dt 



(16) 
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TABLE I: Parameter set for the modulated initial conditions. The column "Parameter set" gives the names of the parameter 
sets and the figures for which a parameter set was used. 



Parameter set 


lattice size 


w (range of fluctuations) ei £2 n m 


e (energy density) 


M-F4 (FigfjTa)) 


4 3 


7.1 x 10" 2 - - 


7.15 x 10" 3 


M-B4 (FigOJb)) 


4 3 


9.7 x 1(T 2 1(T 4 3 3 


7.13 x 10~ 3 


M-Bf4-1 (FigOJc) and HI) 


4 3 


1(T 4 9.7 x 1(T 2 1(T 4 3 3 


7.13 x 1(T 3 


M-f4 (FigGLTc')) 


4 3 


10" 4 - - 


1.379 x 1(T 8 


M-Bf4-2 (Fig© 


4 3 


10" 4 6.859 x 10~ 2 10" 4 3 3 


3.529 x 1(T 3 


M-Bffi (Fig© 


8 3 


10" 4 9.7 x 10" 2 10" 4 3 3 


3.528 x 1(T 3 











D. Choice of initial conditions 



the gauge field configuration created in the collision be- 



As mentioned in the Introduction, realistic initial con- 
ditions for high-energy heavy-ion collisions are based on 



TABLE II: Parameter set for constant A initial conditions. 



Parameter set lattice size w (width of fluctuations) magnetic field B e (energy density) 



C-Bf4 (Fig|3]and|BJ| 
C-Bf8 (Fig© 



5 x 1CT 4 
5 x f(T 4 



0.15 
0.15 



1.45 x 10" 2 
1.45 x 10~ 2 



tween two color-glass condensates, in which both the 
chromomagnetic and chromoelectric fields are parallel to 
the collision axis. These gauge field configurations serve 
as the starting point of the glasma. 

For "modulated initial conditions" (denoted by "M-" 
in the names of the parameter set), the spatial compo- 
nents of the gauge fields are given by an oscillating back- 
ground field plus small fluctuations 



A1{r) =8., 



[2xTT 

eisin U^ 

. / 2nxn\ . ( 2mzir\ 



-t2 sm 



\ N x J 



sm 



Nz J 



Ef(f) = 0. 



+ Vt(r r ), (17) 
(18) 



Here rjf (r) is a random number that has indices for color 
a and spatial direction i (i = 1, 2, 3). 

r/f (f) is a uniformly distributed random number rang- 
ing from —to to to, where to is the amplitude of the 
noise. The superscript a denotes an adjoint color index 
(a = 1, . . . , (iVf — 1)). In the present study, we study 
SU(iV c = 2). N X (N Z ) is the lattice size in x(z) direction. 
Recall that we use the temporal gauge and hence Aq = 0. 
The parameters for the modulated initial condition are 
summarized in Table HI We also include initial conditions 
with fluctuations only (i.e. without background field), 
namely M-F4 and M-f4, for comparison. 

When T)f(r) = 0, the initial magnetic fields are written 
as 



BZ 



Bt 



d y A a z 



d z A a y 



gi ahc A\A\ 



-d z A% 



— 2m— — e2 sin ■ 
d x A% - d x A1 



2nxiT 



2mz7r 



■cos ■ 



d x A a y 

2n 



d v A%- 

2x-k 



N x 

- ge abc A b y A c z 

■ge abc AlA c z 



N> 



d x A a y 



— ei cos 

N x N x 



2mr 2nxiT . 
69 cos sm 

N x N x 



2mzTT 

N z 



(19) 

(20) 

(21) 
(22) 

(23) 



Therefore, the background magnetic fields have x and z 
components, and are modulated along the x and z direc- 
tion. 

The other initial condition adopted here is the 
"constant-^ initial condition" (denoted by "C-" in the 
names of parameter sets), where constant gauge fields 
are used to produce constant magnetic fields (see |26j). 



where g is the coupling constant. The vector potential 
gives the following magnetic field configuration in SU(2), 



B 1 



-F 1 
yx 



9e lbc A b y A x 



9\A y A x — A y A x ) — B 



(26) 



where the other components of the magnetic field are 
zero. Since g is an irrelevant parameter in CYM, we set 
g = 1 throughout the paper. The parameters for the 
constant- A initial condition are summarized in Table HIl 

The constant-A initial condition gives a constant chro- 
momagnetic field in one direction, and can be regarded 
as the simplest modeling of the initial glasma. The ei- 
dominant case of the modulated initial condition (ei 3> 
£2) also simulates the initial condition of the glasma with 
an additional amplitude oscillation in the x direction: 
The color-magnetic field aligns almost in the z direction, 
and a small longitudinal oscillation of the background 
field is introduced by the small value of €2- These features 
of the modulated initial conditions are advantageous for 
the investigation of Nielsen-Olesen instabilities [12(- In a 
strong color-magnetic field, the modulation in the color 
magnetic field direction is most unstable [1J, [lj] . 

While the color-electric field is absent in the present 
initial conditions, the non-linear coupling between the 
chromo-magnetic fields and fluctuations will generate 
chromo-electric fields very fast. It would, of course, be 
desirable to incorporate chromo-electric fields in the ini- 
tial condition for a more comprehensive analysis of thcr- 
malization. However, this is technically more involved 
and beyond the scope of the present work, because one 
would have make sure that the initial condition satisfies 
Gauss' law for zero color charge. 



III. CHAOTICITY WITH GLASMA-LIKE 
INITIAL CONDITIONS 

A. Time evolution of Dee and Def 

Here, we show the time evolution of the distances 
between two trajectories stemming from two adjacent 
points: The parameter sets used for these analyses are 
summarized in Tables HI and HIl Figures QJa), (b), (c) and 
(c') show Dee and Def as functions of time for several 
different initial conditions, M-F4, M-B4, M-Bf4-1 and M- 
f4 in Table HI respectively. The two (adjacent) starting 
points arc related to each other by 



At(7) 



(8 i 2S a3 + S l3 S a2 )VB/g, 
0, 



(24) 

(25) 



A?(t = 0, f) 
E?(t = 0, r) 



1.005A?(i 
1.005£?(t 



0,r), 

0, r) . 



(27) 
(28) 
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FIG. 1: Time evolution of Dee and Def with the initial conditions of (a) only fluctuations (parameter set M-F4 in Table [TJ> , 
(b) only background field (no fluctuations, M-B4), (c) background plus tiny fluctuations (M-Bf4-1), (c') only tiny fluctuations 
whose amplitude is the same as that of (c) (M-f4). 



We note that a linear increase in the semi-logarithmic 
plots shows an exponential growth of the distance, which 
implies an exponential sensitivity of the trajectories to 
the initial value. This is a typical behavior showing the 
chaoticity of the system, which leads to entropy produc- 
tion when combined with Husimi coarse-graining. The 
initial energy densities in the present cases are approxi- 
mately the same except for the case shown in Fig. He'). 

In Fig. QJa), we show the time-evolution of the dis- 
tances Dee and Dff with the initial condition M-F4, 
where the initial field is given by the fluctuation with no 
background magnetic field. We can see the exponential 
growth of Dee and Def from t = 30, which shows the 
chaoticity of the system. After a time around t = 170, 
these distances become saturated. We find almost the 
same saturation times for Dee and Def- This satura- 
tion property is understood to be due to the fixed total 
energy. 

When the background magnetic fields are present in z- 
direction but the fluctuation is absent in the initial con- 
dition (M-B4), the exponential growth of the distances 



does not manifest itself, and they only show a stationary 
oscillating behavior as shown in FigfTJb). In this case, 
thcrmalization is not expected to occur. 

Figure [IJc) shows the time-evolution of Dee and Dff 
with the background magnetic field plus very tiny fluctu- 
ations in the initial condition (M-Bf4-1). These Z?'s show 
an oscillatory behavior in the very early stage, then they 
start to grow exponentially at a certain time. In this 
setup, the onset time of the exponential growth is about 
£ = 50. 

These results suggest that the solution of CYM solely 
with glasma-likc background field at initial time is un- 
stable, but the instability is triggered only when tiny but 
finite fluctuations are imposed on top of the background 
field at initial time. Figure HJc') is the numerical re- 
sult solely with the initial fluctuations whose amplitude 
is the same as that in Fig. [TJc). We can see that the 
distances do not show increasing behavior at least until 
t = 300. Addition of such a tiny fluctuation is essen- 
tial for the exponential growth when it couples with the 
glasma-likc background field in the initial time, as seen in 
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constant on a time scale similar to that observed for the 
modulated initial condition (see Fig. 1(c)). 

The oscillatory behavior in the initial stage can be un- 
derstood by a linearized analysis of the CYM system, 
which shows that the leading-order time dependence is 
given by a Jacobi elliptic function [26| ■ 

The increase of the distances for these different ini- 
tial conditions indicates that chaotic behavior occurs ir- 
respective of the details of the chosen initial conditions, 
as long as they have some (tiny) random fluctuations on 
some coherent background field. 



B. Time evolution of the Lyapunov spectrum 



FIG. 2: Dependence of Dee and Dee to the various size of 
amplitude of fluctuations, w, from modulated initial condition 
(M-Bf4-1). 
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FIG. 3: Time evolution of Dee and Def from the initial 
condition C-Bf4 where the initial background color-magnetic 
field is constant (C-Bf4). 



Fig. QJc). This shows that the exponential growth of the 
distance between two CYM solutions seen in Fig. QJc) is 
not caused by the fluctuations themselves, but it is an 
inherent instability of the background field that is only 
triggered by the fluctuations. 

We show the dependence of Dff on the ratio of the 
fluctuations to the glasma-likc background field in Fig. [2] 
The larger the relative strength of the fluctuations, the 
earlier the onset time of the exponential growth and the 
saturation time of Dff- This result is of phenomenolog- 
ical importance, because we can predict the thermaliza- 
tion time once the ratio of fluctuation to the background 
field is known. 

Figure [3] shows Dff and Dee for a constant A ini- 
tial condition (C-Bf4). We see that both the distances 
grow with more pronounced oscillations of a larger am- 
plitude than in the initial stage for the modulated initial 
condition, but then again become saturated and almost 



As mentioned before, the KS entropy is defined as the 
sum of all positive Lyapunov exponents and gives the en- 
tropy production rate. Therefore, the spectrum of Lya- 
punov exponents in the late stage is of paramount im- 
portance for entropy production in a CYM system. We 
further expect that the sum of positive intermediate Lya- 
punov exponents (ILEs) corresponds to the entropy pro- 
duction rate at a certain time. If this is the case, the time 
evolution of the spectrum of ILEs describes the time evo- 
lution of the thermalization process. Below we report our 
numerical results for the time evolution of the spectrum 
of positive ILEs. 

Modulated initial condition 

Figure 2] shows the time evolution of the spectrum of 
ILEs for the parameter set M-Bf4-1 in Table Q] The ver- 
tical axis is the magnitude of the positive real ILEs and 
the horizontal axis is the label of the ILEs in descending 
order. Imaginary ILEs are plotted as zeros in these fig- 
ures. The time slices for which the different spectra were 
obtained are shown in the upper right-hand corner of the 
figures. 

The figures show that the spectrum has a step-like 
shape at t = 0, which implies the existence of degen- 
erated, unstable, exponentially growing modes from the 
very beginning of the evolution. More importantly, the 
number of the positive Lyapunov exponents is a finite 
fraction of the total number of degrees of freedom. This 
behavior actually could be expected, because the ILE's at 
t = coincide with the local Lyapunov exponents (LLEs) 
defined as the eigenvalues of the Hessian, and thus the 
existence of positive ILE's at t = should reflect that 
of the unstable modes revealed in linear response theory 
0. 

The number of positive ILEs or unstable modes on the 
4 3 lattice is around 60 at t — and increases with time, 
while the magnitudes of the largest ILEs decrease dur- 
ing the time period of < t < 100. From a comparison 
with Fig.QJc), we find that the exponential growth of the 
Z?'s starts well after the spread of the spectrum of unsta- 
ble modes. This observation suggests that a nonlinear 
analysis including mode-mode coupling is necessary for 
the investigation of entropy production. The spectrum 



0.5 
0.45 

0.4 
0.35 

0.3 
0.25 

0.2 
0.15 

0.1 

0.05 





t=0 + 
t=40 



M-Bf4-1 

t (energy density: 0.00713) 



t=70 + 
t-100 

t=150 * 



M-Bf4-1 

h^ (energy density: 0.00713) 




FIG. 4: Time evolution of ILE for the "modulated initial condition" with M-Bf4-1 in Table U(V = 4 3 
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FIG. 5: Time evolution of ILE for the "modulated initial condition" with M-Bf4-2 in Table [J(V = 4 3 



at saturation contains a finite number of positive ILEs. 
In fact, the number of positive ILE's now exceeds 200, 
which constitutes a macroscopic number, comparable to 
the total number of modes, 1152. We conclude that the 
KS entropy is definitely nonzero and that entropy pro- 
duction is a sustained property of the glasma in CYM. 

Figure [5] shows the time evolution of the spectrum of 
ILEs for the parameter set M-Bf4-2, where the initial 
energy density is about half of that for M-Bf4-1. The 
feature of the spectrum is the same as that for M-Bf4- 
1: the step- like shape at the initial time and the time 
dependence of the ILE spectrum. A difference is that the 
magnitude is smaller than for M-Bf4-1 at any point along 
the horizontal axis. This fact shows that KS entropy is 
larger for a larger energy density. We will show the value 
of the KS entropy in Table Hill 

Figure [5] shows the evolution of the ILE spectrum in a 
larger (8 3 ) volume, with the parameter set M-Bf8 shown 
in Table HI The plot style is the same as in Fig. |4] The 
spectrum again starts from a somewhat discontinuous 
and step-like shape at t = and then becomes a smooth 
shape quite similar to that found on the smaller (4 3 ) lat- 
tice. The transition between the two regimes occurs at 
t ~ 50. It should be noted that the total number of de- 
grees of freedom, ISA'' 3 for a N 3 lattice, is quite different 



for the two lattices. Again, the Lyapunov spectrum be- 
comes stable after t ~ 100, containing a large number 
of positive ILE's, approximately 4500, which constitutes 
a macroscopic fraction of the total number of degrees of 
freedom, which is 9216. The similarity with the results 
for the smaller lattice extends to quantitative details. For 
example, the value of the largest ILE is near 0.08 for both 
lattices, and the values of the ILEs for the 17% most un- 
stable modes exceed 0.01 in both cases. 

In FigEl the ILE spectra for M-Bf4-2 and M-Bf8 are 
compared at late time, t = 150, where the spectrum is 
approximately stable. We can see that the magnitude of 
the ILEs for the 8 3 lattice is always larger than for 4 3 
at any point along the horizontal axis. In addition, the 
spectra coincide for large values on the horizontal axis. 
Based on this analysis of volume-dependence, we expect 
that the ILE spectrum stays finite in infinite volume. 

In Table IIII1 the time dependence of the KS entropy 
divided by volume, s KS = S K s/V = Sa>o^/^' is 
shown for each initial condition, M-Bf4-1, M-Bf4-2, and 
M-Bf8. Note that M-Bf4-2 and M-Bf8 have the same 
initial energy density. For the initial condition M-Bf4-1, 
the value of sks = Sks/V is about 0.23 at the initial 
time, drops slowly with time and finally settles around 
0.12 at t ~ 100. Thus, the KS entropy is surely positive 
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FIG. 6: Time evolution of ILE for the same initial condition as in Fig. [5] except for the lattice size (8 in this figure). The 
parameter set used in the calculation is M-Bf8 in Table [I] 
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FIG. 7: Dependence of ILE spectrum on lattice size. The 
spectrum is plotted against the label of the ILEs divided by 
the total number of modes in descending order of ILE. The 
red-solid line shows the spectrum on the 4 3 lattice with the 
initial condition M-Bf4-2, and the green-dotted line shows 
that on the 8 3 lattice with the initial condition M-Bf8. These 
two initial conditions have the same energy density. 



definite even at late times, implying that the entropy is 
produced. A similar behavior is obtained for the initial 
condition M-Bf4-2 which has an energy density about 
half of that of M-Bf4-1; sks starts from a value around 
0.15, and the stable final value around 0.08 is reached at 
t ~ 150. The larger the energy density, the larger the 
value of the KS entropy. In the previous study, Ref. [25| , 
the authors found that sks scales as the fourth root of 
energy density, e 1 ' 4 , for random initial magnetic fields. 
For the initial conditions M-Bf4-1 and M-Bf4-2, the ra- 
tio of KS entropy, 0.122/0.079 ~ 1.54, is not so close 
to the fourth root of the ratio of the energy density e, 
(0.00706/0.00353) 1 / 4 ~ 1.19. As the setting in this case 
is not isotropic, it is not clear that all dimensioned pa- 
rameters should simply scale with powers of e 1 ' 4 . It could 
also be that the volume studied here (4 3 ) is insufficient 
for a reliable scaling analysis. 



For the larger volume, M-Bf8, sks is about 0.14 at 
initial time and about 0.10 at t = 150. sks surely sur- 
vives in a larger volume, and entropy production in the 
infinite-volume limit can be expected. 

Constant-A initial condition 



We show the time evolution of the ILE spectra ob- 
tained from the constant- A initial conditions with tiny 
fluctuation on 4 3 and 8 3 lattices in Figs. [8] and H3 respec- 
tively. At t = 0, the spectra have the step- like structure 
similar to that from the modulated initial condition, ex- 
cept for the detailed structures. At later time, this dis- 
continuity of the structure disappears and the spectrum 
becomes smooth. After about t = 100, the spectrum be- 
comes stable, and again the fraction of positive ILEs is 
finite. The behavior on the 8 3 lattice is quantitatively 
quite similar to that on the 4 3 lattice, when the horizon- 
tal axis is rescaled by the total number of modes. 

Comparison of the two cases we have studied, modu- 
lated and constant- A initial conditions, which are chosen 
to mimic the glasma initial condition, reveals that the 
ILE spectra at late times are quite similar: Except for 
a few modes that remain very unstable (A > 0.1) in the 
case of the constant- A initial condition, the largest ILE's 
take values around 0.08, and the modes at the 17 th per- 
centile of all modes (the 200 th and 1600 th modes for the 
4 3 and 8 3 lattices, respectively) have ILEs around 0.01. 
In both cases also a much larger fraction of all modes 
becomes unstable at late times than in the initial phase 
of the evolution. The agreement suggests that entropy 
production occurs in the CYM dynamics for any glasma- 
like initial condition and that a large part of the entropy 
production is caused by the chaoticity of the gauge field 
dynamics manifested after some time, rather than by the 
instability of the specific gauge field configuration in the 
initial stage. 
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FIG. 8: Time evolution of the intermediate Lyapunov exponent spectrum for the "constant A initial condition" on the V — 4 3 
lattice. The parameter set used in the calculation is C-Bf4 from Table iHl 
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FIG. 9: Time evolution of intermediate Lyapunov exponent spectrum for the "constant A initial condition" on the V 
lattice. The parameter set used in the calculation is C-Bf8 from Table llll 
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TABLE III: Time dependence of Kolmogorov-Sinai' entropy 
divided by volume for initial conditions, M-Bf4-1, M-Bf4-2 
and M-Bf8. e denotes the energy density. 



IV. SUMMARY AND CONCLUDING 
REMARKS 



We have studied the thermalization process of classi- 
cal Yang-Mills fields starting from semi-realistic glasma- 
like initial conditions characteristic for relativistic heavy- 
ion collisions. We focused on the chaotic behavior and 



entropy production of the classical field theory. For 
this purpose we determined all Lyapunov exponents for 
the theory discretized on a cubic spatial point lattice. 
The sum of all positive Lyapunov exponents, i.e. the 
Kolmogorov-Sinai entropy, gives the rate of entropy pro- 
duction. 

We considered two types of glasma-like initial condi- 
tions. One is the "modulated initial condition" , where 
initial chromomagnetic fields Bf are spatially modulated 
along the z and x axes. In the absence of fluctuations, 
this field configuration is color- independent. The other 
initial condition is the "constant- A initial condition" , 
where the gauge potentials Af and chromo-magnetic 
fields B? are both constant, but color-dependent, and 
only the ^-component of the chromomagnetic field in one 
color direction is nonzero. For both types of initial con- 
ditions we added small random fluctuations to describe 
the noise present in the glasma fields due to the quantum 
fluctuations in the two colliding nuclei. 

We have found that when the gauge field is given by 
the modulated initial condition without any fluctuations, 
the distance between two trajectories starting from ad- 
jacent points shows oscillatory behavior and exponential 



11 



growth is absent. When small fluctuations are present 
on top of the background field, the initial oscillatory be- 
havior of the distance terminates after a short time and 
the exponential growth of the distance becomes robustly 
visible. The onset time of exponential growth depends 
on the ratio of the amplitude of the fluctuation to the 
strength of the chromomagnetic field. The fact that we 
do not see any exponential growth in the absence of initial 
fluctuations demonstrates that our numeric treatment in- 
troduces no fluctuations of relevant size. 

For the modulated initial condition with small fluctu- 
ations positive Lyapunov exponents are present from the 
beginning. This is a reflection of the unstable modes re- 
vealed in the linear response analysis of glasma fields. We 
found that the number of modes with positive Lyapunov 
exponents increases substantially during time evolution, 
and that the Lyapunov spectrum reaches a stable shape 
when a macroscopic fraction of the modes has positive 
Lyapunov exponents. This implies that a bulk number 
of field modes contributes to the KS entropy and that 
the entropy production rate per unit volume is non-zero 
in CYM. 

The dependence of the Lyapunov spectrum on the size 
of the spatial volume is found to be weak if the mode 
number is rescaled by the total number of modes, which 
suggests that extensive entropy production occurs in the 
infinite volume limit in CYM as first found by Bolte et 
al. (27J|. 

The chaotic behavior of CYM with the noisy 
"constant- A initial condition" is similar to that observed 
for the modulated initial condition with some modifica- 
tions during the initial phase of the time evolution. After 
an initial increase of the distance accompanied with some 
oscillatory behavior, the distance shows an exponential 
growth, and the time evolution of the spectrum of ILEs 
shows a qualitatively similar behavior as that seen for 
the modulated initial condition. The number of positive 
Lyapunov exponents is again found to be large, corre- 
sponding to many unstable bulk modes, and the entropy 
production rate grows approximately linear with volume 
also for the constant A initial condition. 

The important conclusion to be drawn from our 
present analysis is that entropy production occurs in a 
robust way at the classical field level starting from semi- 
realistic initial conditions. While the classical field de- 
scription of the glasma becomes invalid when thermal 
equilibrium is approached, the pre-equilibrium dynamics 
of the glasma can be simulated by classical gauge field 
equations, which allows to estimate the thermalization 
time. We therefore focused our study on the intermediate 
time regime, where the ILEs characterize the dynamical 
instability of the gauge field. Our study suggests that 
the magnitude of the delay time before the start of lin- 
ear entropy growth depends substantially on the chosen 
initial condition while the entropy growth rate itself is 
affected at most mildly. Substantially more systematic 
studies are needed to decide whether there exist realistic 
scenarios for which the thermalization time is as small as 



1 fm/c, which is the value advocated by the comparison 
of hydrodynamical simulations with heavy ion data. 

As the gauge field approaches equilibrium and quan- 
tum effects become important by providing a physical 
cut-off to the ultraviolet divergences of the classical ther- 
mal field theory, the evolution of the glasma can be de- 
scribed by viscous hydrodynamics. It is an interesting 
question whether the two effective theories of the dynam- 
ics of an equilibrating quark-gluon plasma can be directly 
joined or whether there is a need for an additional for- 
malism interpolating between these two regimes. 
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Appendix : Proof of the gauge invariance of 
Lyapunov spectra 

We here show that the local/intermediate Lyapunov 
exponents are gauge invariant under the remaining 
gauge transformations in temporal gauge, i.e., time- 
independent gauge transformations. The Hessian Hit) 
in a concrete expression is 



H(t) 




(1) 



where {H AA (t)) iax j by = 5 2 H/5A'?{x,t)SA b j {y 1 t) is a 
second derivative in terms of gauge fields. Gauge 
fields Af are transformed according to A'f(x,t) = 
(p,(x)Ai(x,t)) a + W°-(x,t), with a time-independent or- 
thogonal matrix O(af), i.e., £l(x)n, T (x) = 1. Wf(x, t) is 
a possible term independent of gauge fields. 
Paying attention to the chain rule, 



SAUx,t) 



6A'f(2,t) 5A b {x,t)6A"*(x,t) SA b (x,t) 



[n T {x)) l 



the transformation property of Uaa for a time- 
independent gauge transformation Q(x) in the adjoint 
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representation of SU(N) is found to be 

H AA (t) = 5 2 H/5A1{x,t)5A){y,t) 
^H' AA (t) = 5H 2 /5A la l {x 1 t)5A'){y,t) 

= (0,H AA (t)Q )iax,jby (2) 

(Note that the Hamiltonian H is gauge invariant.) Defin- 
ing an orthogonal martix 



VL{x) 



n(x) o 



(3) 



the gauge transformation of %{t) can be expressed as 
H(t)^H'(t)=TiH{t)TL T 



(4) 



in a similar manner. Taking orthogonality and time- 
independence of Q(x) into account, the local and inter- 
mediate Lyapunov spectra of the Hessian rt{t) are easily 



shown to be gauge invariant. In fact, a time-ordered 
product of the Hessian 



U(t,t + At) =T 



transforms as 



/ rt+At 

exp / H(t + t')dt' 



(5) 



U(t,t + At) -+ U'(t,t + At) 



T 



exp 



nr 




t+At 



nn{t + t')n dt' 



t+At 

I H{t + t')dt' 



— T 





nu(t, t + At)n 



T 



(6) 



and gives the same (gauge invariant) intermediate Lya- 
punov spectra. 
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